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1.  Introduction 

The  objective  of  this  research  effort  is  to  develop  and  test  a  methodology  for 
regional  seismic  event  identification  which  utilizes  2-dimensional  signal  processing 
and  treats  the  seismogram  as  an  image.  The  methodology  seeks  to  exploit  the 
information  content  of  an  entire  regional  seismogram,  and  is  particularly 
appropriate  for  small  events,  multiple  events,  and  events  with  blocked,  attenuated, 
or  missing  phases. 

The  motivation  for  this  research  comes  from  our  observation  that  the  spectral 
features  which  distinguish  quarry  explosions  from  naturally  occurring  earthquakes 
can  be  seen  throughout  the  entire  seismogram.  These  features,  which  include  the 
spectral  complexity  and  spectral  energy  distribution  (P  vs.  S),  are  readily  apparent  in 
the  time-frequency  (TF)  representation  of  the  seismogram1.  The  TF  distribution  is 
particularly  helpful  in  the  analysis  of  events  which  have  highly  attenuated,  blocked, 
or  missing  phases. 

A  number  of  investigators  have  pointed  out  the  utility  of  the  TF  distribution 
for  seismic  analysis  (Smith,  1989;  Hedlin  et  ah,  1989).  However,  TF  plots  have  not 
been  widely  used  or  accepted  in  the  seismic  community.  There  are  a  number  of 
reasons  for  this  situation:  the  difficulty  of  manipulating  and  operating  on  2-D  images 
us.  1-D  time  series  or  spectra;  the  lack  of  quantifiable  measures  of  the  features  seen 
on  TF  plots;  the  lack  of  simple  relations  between  TF  images  and  the  physics  of  the 
source  and  propagation  effects;  and  the  lack  of  a  simple  image  discriminator  which 
would  allow  automated  processing  of  TF  images. 

During  the  course  of  our  research  on  the  application  of  artificial  neural 
networks  (ANN’s)  to  seismic  event  identification,  we  have  employed  TF  plots  for 

1  The  time-frequency  distribution  goes  by  many  names.  Examples  include  spectrogram  and  sonogram. 
We  prefer  the  general  term  time-frequency  distribution,  because  it  includes  many  different  types  of 
distributions,  such  as  Wigner-Ville,  and  special  window  functions  such  as  the  cone  kernel.  A 
spectrogram  or  sonogram  implies  a  specific  type  of  distribution. 
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the  purposes  of  data  review  and  phase  timing.  This  use  has  led  us  to  be  able  to 
perform  simple  visual  classification  of  TF  images  which  closely  match  the 
classification  capability  of  our  1-D  ANN-based  processor.  However,  the  TF  image 
appears  to  excel  in  situations  of  low  SNR,  mixed  events,  or  events  with  highly 
attenuated  phases.  Preliminary  research  has  shown  that  we  can  exploit  the 
information  content  of  a  TF  distribution  by  first  transforming  the  TF  image  into  the 
wavenumber  domain  via  a  2-D  FFT,  and  then  use  a  combination  of  image 
processing  and  simple  backpropagation  ANN's  to  extract  a  number  of  features 
which  completely  describe  the  spatio-temporal  nature  of  the  source.  These  features 
can  then  be  input  to  any  number  of  classifiers  to  perform  the  final  event 
identification. 

This  report  documents  our  study  to  test  the  utility  of  these  concepts  for 
regional  seismic  event  identification.  The  study  focuses  on  events  from  the 
NORESS,  Vogtland,  and  Lubin  ground  truth  databases  used  in  our  previous  studies 
of  regional  discrimination  (Pulli  and  Dysart,  1993). 
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2.  Signal  Processing  and  Parameter  Extraction 


2.1  Step  1:  Time-Frequency  Transformation 

Step  1  of  our  analysis  involves  the  transformation  of  the  original  time  series 
data  (single  station  or  array)  into  the  time-frequency  (TF)  domain.  There  are 
numerous  schemes  for  implementing  the  TF  transformation  (Cohen,  1989),  the 
simplest  of  which  is  to  slide  a  tapered  window  down  the  time  series  and  compute 
successive  FFT's.  The  result  can  be  displayed  in  a  number  of  ways,  but  we  prefer  a  2- 
dimensional  color-coded  representation  for  simplicity.  Figure  1  shows  the  TF  images 
of  two  demonstration  events  recorded  at  the  NORESS  array.  The  event  at  the  top  is  a 
quarry  blast  at  the  Titania  Mine,  southwest  of  NORESS.  Note  the  spectral  banding 
seen  throughout  the  TF  image.  The  event  at  the  bottom  is  a  nearby  earthquake.  No 
spectral  banding  is  present  in  its  TF  image.  The  TF  image  of  the  earthquake  is 
dominated  by  energy  around  the  Lg  wave. 

2.2  Step  2:  Image  Equalization 

Image  equalization  refers  to  the  process  by  which  the  images  are  enhanced, 
corrected,  or  detrended  to  bring  out  specific  features  of  the  image.  An  analogy  in  1- 
dimensional  seismic  signal  processing  would  be  the  correction  of  seismic  spectra  for 
/“2  spectral  decay  or  attenuation.  Three  different  image  equalization  schemes  were 
applied  to  the  TF  distributions. 

Equalization  1  treats  the  TF  image  as  a  sequence  of  power  time-series 
(essentially  narrow  bandpass  filtered  seismograms,  one  frequency  at  a  time)  and 
applies  noise-spectral  equalization  to  the  rows  (time  series)  of  the  image.  This 
equalization  emphasizes  the  relative  energy  in  each  frequency  bin  as  a  function  of 
time,  and  tends  to  enhance  the  Pn/Lg  spectral  ratio  information  in  the  whole 
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seismogram.  It  also  enhances  the  spectral  banding  seen  in  the  ripple-fired  quarry 
blasts  (Smith,  1989;  Hedlin  et  ah,  1989). 

Equalization  2  treats  the  TF  image  as  a  sequence  of  power  spectra  and  applies 
noise-spectral  equalization  to  the  columns  (spectra)  of  the  image  to  obtain  deviations 
from  the  trend  of  the  power  spectra.  This  equalization  emphasizes  the  spectral 
banding  seen  in  the  ripple-fired  quarry  blasts  and  de-emphasizes  the  Pn/Lg  spectral 
ratio  information. 

Equalization  3  is  similar  to  equalization  1  in  that  it  operates  on  the  rows  (time- 
series)  of  the  TF  distribution,  but  it  this  case  it  determines  deviations  in  the  trends  of 
the  time  series.  This  effectively  removes  the  seismic  phase  amplitude  ratio 
information  and  the  frequency  banding,  leaving  whatever  is  left  to  dominate. 
Hence,  this  equalization  is  a  residual  estimate  of  the  TF  iamge. 

2.3  Step  3:  Wavenumber  Transformation 

In  Step  3,  the  TF  image  is  transformed  into  the  wavenumber  domain  using  a 
2-D  FFT.  The  wavenumber  representation  of  the  demonstration  events  is  shown  in 
Figure  2.  A  variety  of  options  exist  at  this  stage,  including  the  use  of  alternative  2-D 
images  such  as  cepstra  vs.  time,  Lg/Pn  spectral  ratios  vs.  time,  binary  images,  and 
higher  dimensional  images. 

2.4  Step  4:  Image  Model  and  Image  Parameters 

A  key  element  of  our  methodology  is  the  use  of  a  parametric  model  to 
characterize  the  power  spectrum  of  an  arbitrary  image.  The  model  we  have  used  is  a 
stochastic  surface  model  originally  developed  by  Goff  and  Jordan  (1989)  to 
characterize  bathymetric  surfaces.  The  power  spectral  form  of  the  model  is  given  by 
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m= 


47wH2 

M 


[uHk)+ lp" 


(1) 


where  H  is  the  RMS  height  (energy),  k  is  the  wavenumber,  Q  is  a  matrix  given  by 


Q  = 


k2  sin2  ( 0)  +  kl  cos2  ( 0) 
{k2  -£2)sin2(0)cos2(0) 


[k2  -&2)sin2(0)cos2(0) 
k2  sin2  ( 6)  +  k2  cos2  ( 6) 


(2) 


and  0  is  the  azimuth.  The  exponential  term  v  is  related  to  a  constant  asymptotic  slope 
(B)  of  the  power  spectrum,  or  fractal  dimension  (D),  each  of  which  implies  a  surface 
whose  amplitudes  scale  predictably  with  wavenumber.  B  and  D  are  related  to  v  by 


B  =  — 2(  t>  + 1) 


(3) 


D-3-v 


These  parameters  quantify  features  such  as  lineations,  spatial  orientation, 
characteristic  lengths,  and  roughness. 

In  summary,  the  model  of  the  wavenumber  image  is  characterized  by  5 
features: 

•  H:  the  RMS  height,  or  energy  of  the  wavenumber  spectrum 

•  0:  the  azimuth  of  symmetry 

•  ks:  the  characteristic  wavenumber  parallel  to  the  axis  of  symmetry 

•  kn:  the  characteristic  wavenumber  normal  to  the  axis  of  symmetry 

•  B:  the  asymptotic  slope 
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In  addition  to  these  5  model  parameters,  we  extracted  an  additional  19  image 
parameters  using  techniques  from  image  processing  (such  as  the  Karhunen-Loeve 
transform).  This  transform  determines  the  covariance  structure  and  higher-order 
moments  of  the  image  by  assuming  that  the  image  represents  a  pseudo-probability 
density  function.  The  covariance  of  the  zero-mean  image  is 

COV[P(bc,ky)]  =  {S"  M  W 

[a21  o22  J 


where 


■s„=Z 


( kxi  -  kx)2 
n  —  1 


Y,P(kx,ky) 

j= i 


{kyi  -  ky)2 
n  —  l 


I 


P(kx,ky) 


■5,2=1 


X 


(kxt  ~  kx){ky]  -  ky) 
n  —  l 


- P(kx,ky ) 


(5) 


The  eigenvalue  problem  is  solved  by  the  singular  value  decomposition  (SVD)  of  the 
covariance  matrix 

svd(cov)  =  f/j^  °2|yr  (6) 


Although  not  required  in  the  2X2  case,  the  SVD  method  is  included  by  our  analysis 
to  allow  parameter  extraction  from  higher-dimensional  images. 

The  parameters  determined  by  this  analysis,  as  well  as  their  names  as  used  in 
subsequent  plots,  are  shown  in  the  following  table. 
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Table  1.  List  of  parameters  extracted  from  wavenumber  images. 


Parameter 

Definition 

Parameter  Name  on  Plots 

eil 

1st  Eigenvalue,  i=l 

emajor 

e12 

1st  Eigenvalue,  i=2 

emajor-2 

e13 

1st  Eigenvalue,  i=3 

emajor-3 

e21 

2nd  Eigenvalue,  i=l 

eminor 

e22 

2nd  Eigenvalue,  i=2 

eminor-2 

e23 

2nd  Eigenvalue,  i=3 

eminor-3 

kurti 

Image  kurtosis,  i=l 

curt 

kurt2 

Image  kurtosis,  i=2 

curt-2 

kurt3 

Image  kurtosis,  i=3 

curt-3 

eil/e21 

Eigenvalue  ratio,  i=l 

emajor/eminor 

ei2/e22 

Eigenvalue  ratio,  i=2 

emajor-2/eminor-2 

e13/ e23 

Eigenvalue  ratio,  i=3 

emajor-3/ eminor-3 

en  x  e2l 

Eigenvalue  product,  i=l 

emajor  x  eminor 

e12  x  e22 

Eigenvalue  product,  i=2 

emajor-2  x  eminor-2 

ei3  x  e23 

Eigenvalue  product,  i=3 

emajor-3  x  eminor-3 

4eu+el\ 

Eigenvalue  norm,  i=l 

eigen_norm 

Ve  12  +  e22 

Eigenvalue  norm,  i=l 

eigen2_norm 

Vg13  g23 

Eigenvalue  norm,  i=l 

eigen3_norm 

e13/ e22 

Eigenvalue  ratio,  i-3,  j=2 

emaj  or-3  /  emaj  or-2 

Kurt3  /  Kurt2 

Kurtosis  ratio 

curt-3 /curt-2 

ks 

Characteristic  length 

1/ws 

kn 

Characteristic  length 

1/wn 

H 

RMS  height 

rms 

B 

Spectral  slope 

hexp 

Frequency  (Hz)  ,  Nanometers 


-10  0  10 


Wavenumber 


-10  0  10 


Wavenumber 


Figure  2.  Wavenumber  spectra  of  Titania  mine  blast  (top)  and  coastal  earthquake  (bottom) 
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3.  Results  for  NORESS,  Vogt  land,  and  Lubin 


3.1  NORESS 

The  methodology  just  described  was  first  applied  to  83  events  from  our 
NORESS  ground  truth  database  described  in  previous  documents  (e.g.  Pulli  and 
Dysart,  1994).  Of  these  83  events,  35  were  earthquakes  and  48  were  explosions.  The  24 
image  parameters  were  extracted  for  all  of  the  events  using  each  of  the  three  image 
equalization  schemes.  Values  of  these  parameters  are  shown  graphically  in  Figures 
3-8. 

It  is  clear  from  these  plots  that  some  of  the  image  parameters  are  quite  good  at 
separating  event  types,  such  as  curtosis  for  equalizations  1  and  2.  However,  many 
parameters  provide  little  or  no  capability  for  separating  event  classes.  An  effective 
means  of  condensing  and  evaluating  these  parameters  is  the  calculation  of  the 
Mahalanobis  Distance,  which  is 


(fli  -/*2)2 

2  .  _2 
CT,  +( J2 


(7) 


Here,  |i  and  cr  are  the  means  and  variances  respectively  of  the  parameter  values  for 
each  class  of  events  (1  and  2).  Parameters  which  provide  good  identification 
capability  have  large  differences  in  mean  values  and  small  variances,  and  hence 
large  Dp's.  The  best  parameters  from  equalization  1  are: 

•  20:  eminor  3 

•  1:  curtosis 

•  10:  emajor-2  /  eminor-2 

•  18:  eminor 
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For  equalization  2,  the  best  parameters  are: 


•  2:  curtosis  2 

•  1:  curtosis 

•  4:  curtosis-3  /  curtosis-2 

•  22:  rms 

•  19:  eminor  2 

For  equalization  3,  the  best  parameters  are: 

•  15:  emajor-3  x  eminor-3 

•  19:  eminor  2 

•  4:  curtosis-3  /  curtosis-2 

•  22:  rms 

•  2:  curtosis-2 

3.2  Vogtland 

The  Vogtland  dataset  (see  Grant  et  ah,  1994;  Pulli  and  Dysart,  1994)  consists  of 
23  events,  9  earthquakes  and  14  explosions.  Parameter  values  for  the  three  image 
equalization  schemes  are  shown  in  Figures  12-17.  Mahalanobis  Distances  are  plotted 
in  Figures  18-20.  The  best  parameters  from  equalization  1  are: 

•  18:  eminor 

•  20:  eminor-3 

•  1:  curtosis 
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•  15:  emajor  -3  x  eminor-3 

•  4:  curtosis-3  /  curtosis-2 

•  22:  rms 

For  equalization  2,  the  best  parameters  are: 

•  18:  eminor 

•  1:  curtosis 

•  23:  wn 

•  24:  ws 

•  21:  hexp 

For  equalization  3,  the  best  parameters  are: 

•  9:  emajor  2 

•  5:  eigen2_norm 

•  13:  emajor-3  /  emajor-2 

•  11:  emajor-2  x  eminor-2 

3.3  Lubin 

The  Lubin  dataset  consists  of  30  "mining-induced  earthquakes”.  Image 
parameters  for  the  three  equalizations  are  shown  in  Figures  21-26.  Since  there  is 
only  one  source  type  in  this  dataset,  Mahalahobis  Distance  calculations  are  not 
appropriate,  but  a  comparison  of  the  image  parameters  with  those  of  the  nearby 
Vogtland  events  is  of  interest.  We  will  take  the  best  parameters  from  equalization 
and  note  under  which  class  the  Lubin  events  fall. 
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The  best  parameters  from  equalization  1  are: 


•  18:  eminor  -  Lubin  events  fall  primarily  under  explosion  class 

•  20:  eminor-3  -  Lubin  events  fall  under  both  classes 

•  1:  curtosis  -  Lubin  events  fall  under  earthquake  class 

•  15:  emajor-3  x  eminor-3  -  Lubin  events  fall  primarily  under 
explosion  class 

•  4:  curtosis-3  /  curtosis-2  -  Lubin  events  fall  primarily  under 
explosion  class 

•  22:  rms  -  Lubin  events  fall  under  earthquake  class 
For  equalization  2,  the  best  parameters  are: 

•  18:  eminor  -  Lubin  events  fall  primarily  under  explosion  class 

•  1:  curtosis  -  Lubin  events  fall  primarily  under  explosion  class 

•  23:  wn  -  Lubin  events  fall  primarily  under  explosion  class 

•  24:  ws  -  Lubin  events  fall  primarily  under  explosion  class 

•  21:  hexp  -  Lubin  events  fall  primarily  under  explosion  class 


For  equalization  3,  the  best  parameters  are: 

•  9:  emajor  2  -  Lubin  events  fall  under  both  classes 

•  5:  eigen2_norm  -  Lubin  events  fall  primarily  under  explosion  class 

•  13:  emajor-3  /  emajor-2  -  Lubin  events  fall  under  both  classes 

•  11:  emajor-2  x  eminor-2  -  Lubin  events  fall  primarily  under 
explosion  class 
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Figure  4.  Values  of  image  parameters  13-24  for  the  NORESS  dataset,  equalization  number  1. 
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Figure  6.  Values  of  image  parameters  13-24  for  the  NORESS  dataset,  equalization  number  2. 
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Figure  7.  Values  of  image  parameters  1-12  for  the  NORESS  dataset,  equalization  number  3. 
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Figure  10.  Mahalanobis  distances  for  image  parameters  for  NORESS  database,  equalization  2. 
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Figure  11.  Mahalanobis  distances  for  image  parameters  for  NORESS  database,  equalization  3, 
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Figure  12.  Values  of  image  parameters  1-12  for  the  Vogtland  dataset,  equalization  number  1. 
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Figure  13.  Values  of  image  parameters  13-24  for  the  Vogtland  dataset,  equalization  number  1. 
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Figure  14.  Values  of  image  parameters  1-12  for  the  Vogtland  dataset,  equalization  number  2. 


-25- 


fclW 


fcjW 


fclM 


FsTrij 


0,0  O  > 


^>^<xx 


o°  o  xx  Xn 


%  <p*x 

^3  X  X 


Event  No. 

Event  No. 

Event  No. 

emajor/eminor 

emajorxeminor 

eminor 

%  cP  X 


xx^x 


°  xx 


^  Oflfx  X«V 


°  xxx 
'  CQ3  CE^<  ^xX 
%  °  *  * 


°D  <B  X 

o  o  x 

cP  X  xxx 

O  x 


rQ  A  ' 

Vo  X  X 


x  x  X 

X*X  x 

x  xX<  x 


Event  No. 

Event  No. 

Event  No. 

eminor2 

eminor3 

hexp 

<*%P  X 


->*x  x 


%  o  X 


50  °  X 

00  °  X  Xx 

x  ^x>^x 

X 

5" 

Event  No. 


'JJ  o 

°  o° 
50  %  xv 


0  10  20 
Event  No. 


Figure  15.  Values  of  image  parameters  13-24  for  the  Vogtland  dataset,  equalization  number  2. 
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Figure  16.  Values  of  image  parameters  1-12  for  the  Vogtland  dataset,  equalization  number  3. 
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Figure  17.  Values  of  image  parameters  13-24  for  the  Vogtland  dataset,  equalization  number  3. 
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Figure  18.  Mahalanobis  distances  for  image  parameters  for  Vogtland  database,  equalization  1. 
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Figure  19.  Mahalanobis  distances  for  image  parameters  for  Vogtland  database,  equalization  2. 
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Figure  20.  Mahalanobis  distances  for  image  parameters  for  Vogtland  database,  equalization  3. 
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Figure  21.  Values  of  image  parameters  1-12  for  the  Lubin  dataset,  equalization  number  1. 
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Figure  22.  Values  of  image  parameters  12-24  for  the  Lubin  dataset,  equalization  number  1. 
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Figure  23.  Values  of  image  parameters  1-12  for  the  Lubin  dataset,  equalization  number  2. 
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Figure  24.  Values  of  image  parameters  12-24  for  the  Lubin  dataset,  equalization  number  2. 
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Figure  25.  Values  of  image  parameters  1-12  for  the  Lubin  dataset,  equalization  number  3. 
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Figure  26.  Values  of  image  parameters  12-24  for  the  Lubin  dataset,  equalization  number  3. 
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4.  Summary  and  Conclusions 


The  goal  of  this  research  was  to  determine  if  parameters  extracted  from  the 
time-frequency  (TF)  representations  of  regional  seismograms  could  be  used  to 
identify  seismic  source  types.  To  accomplish  this,  we  extracted  24  image  parameters 
from  the  wavenumber  spectra  of  previously  identified  seismic  events  from  the 
NORESS,  Vogtland,  and  Lubin  ground-truth  databases.  A  total  of  135  events  were 
analyzed. 

A  number  of  the  image  parameters  were  found  to  provide  excellent  source 
identification  capability.  For  example,  the  eigenvalues  of  the  images  and  the  image 
kurtosis  were  found  to  be  useful  for  the  NORESS  and  Vogtland  datasets.  The  Lubin 
events,  which  are  all  "mining-induced  earthquakes",  fall  within  either  the 
earthquake  or  explosion  classes,  which  corroborates  the  results  of  our  1-dimensional 
study  (Pulli  and  Dysart,  1994).  In  general,  we  find  that  the  image  parameters 
extracted  from  the  wavenumber  spectral  images  provide  additional  useful 
information  for  source  identification. 

If  there  is  a  problem  with  the  methodology,  it  lies  in  the  physical 
interpretation  of  the  image  parameters.  Some  parameters,  such  as  the  characteristic 
wavelengths,  reflect  spectral  banding  in  the  TF  distributions.  However,  some  useful 
parameters,  such  as  the  image  eigenvalues,  are  not  readily  related  to  the  source 
physics.  To  overcome  this  difficulty,  a  model  study  should  be  undertaken  so  that 
image  parameters  can  be  extracted  from  TF  plots  of  model  events  with  known 
spectral  characteristics. 
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